Van der Waals interactions in the ground state of Mg(BH 4 ) 2 from density functional 

theory 
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In order to resolve an outstanding discrepancy between experiment and theory regarding the 
ground-state structure of Mg(BH4)2, we examine the importance of long-range dispersive interac- 
tions on the compound's thermodynamic stability. Careful treatment of the correlation effects within 
a recently developed nonlocal van der Waals density functional (vdW-DF) leads to a good agreement 
with experiment, favoring the Q-Mg(BH4)2 phase (P6i22) and a closely related Mn(BH4)2-prototype 
phase (P3il2) over a large set of polymorphs at low temperatures. Our study demonstrates the need 
to go beyond (semi)local density functional approximations for a reliable description of crystalline 
high-valent metal borohydrides. 

PACS numbers: 88.30.rd 63.20.-c 71.20.-b 78.30.-j 
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I. INTRODUCTION 

Light-weight metal borohydrides remain in the spot- 
light of hydrogen-storage research due to their high gravi- 
metric hydrogen content and abundance of the con- 
stituent elements*^ Li, Mg, and Ca tetraborohydridcs 
(M(BH 4 ) n ) hold over 12 wt% of H, but their excessive 
thermodynamic stability presently renders these mate- 
rials unsuitable for practical reversible hydrogen-storage 
solutions^ In an effort to adjust the formation enthalpies 
and the decomposition temperatures, several groups have 
suggested destabilization routes via reactions with other 
hydrides (see e.g. Refs.2r— ) or synthesis of mixed metal 
tetraborohydrides (e.g. Li-Cu) e . Detailed knowledge of 
materials' ground states becomes essential, as the possi- 
bility of new compound formation depends on small free 
energy differences for the phases involved. 

Metal tetraborohydrides are held together primarily 
by a combination of covalent (B-H) and ionic ([M] nS+ - 
[BH4] ) interactions. Packing of the fairly rigid [BH 4 ] 5 ~ 
units and [M] n<5+ ions into crystalline structures de- 
pends strongly on the valency and size of the metal 
atom, with exceptionally complex configurations occur- 
ring for the medium-sized divalent magnesium. Based 
on the experimental data, the ground states below and 
above T — 453 K, namely a and (3, have been assigned 
space groups P6i22 (330 atoms/u.c.} 7 * 8 and Fddd (704 
atoms/u.c.}2, respectively. Strikingly, numerous density 
functional theory (DFT) studies^iSr— have converged on 
a completely unrelated F222 structure as the most stable 
low-T polymorph ] 12 i 14 ' 15 Recently, it has been pointed 
outi& that the lack of a proper description of dispersion 
forces in the standard DFT approach may be a factor fa- 
voring low-density F222 structure over the experimental 
a phase. 

In this study we take a key step towards resolving the 
existing discrepancy: we demonstrate that the weak dis- 
persive interactions, important in other ionic (KC1 and 
KBr) 1 ! and covalent- ionic (Mg(OH) 2 and Ca(OH) 2 )±£ 



systems, indeed play a critical role in defining the 
Mg(BH4)2 ground state. We show that the commonly 
used Ceperley-Alder— or Perdew, Burke, and Ernzerhof 
(PBE)— functionals artificially stabilize configurations 
with unusually high and low density, respectively. In- 
clusion of the dispersive contributions via a nonlocal van 
der Waals density functional (vdW-DF) 21 i 22 or as a semi- 
empirical PBE-D* correction 2 ^— changes the relative 
stability of considered structures and favors the layered 
motif 2 " occurring in related a-Mg(BH 4 ) 2 and Mn(BH 4 ) 2 
prototypes. 

The careful re-examination of Mg(BH 4 ) 2 has become 
possible due to a recent surge of attempts to include vdW 
interactions in the DFTi 27 i 28 The truly nonlocal correla- 
tion functional in the vdW-DF approach has shown good 
transferability for a range of systems reaching from sim- 
ple dimers 2 ^ and physisorbed molecules^ to DNA^i and 
drug design^. An efficient FFT formulation of vdW- 
DF— has allowed us to calculate the T — K relative 
stability for Mg(BH 4 ) 2 structures of unprecedented size 
(up to 330 atom/u.c). Calculation of the Gibbs en- 
ergy vibrational contributions at finite temperatures is 
much more computationally demanding. In order to ex- 
amine the relative stability as a function of temperature 
we have employed a semi-empirical method developed by 
Grimm o 24 ' 25 , that includes the long-range contributions 
via damped pairwise fdmp(R)CeR~ 6 terms at a negligi- 
ble cost compared to standard DFT calculations. The 
method has gained popularity providing an improved de- 
scription of molecular systems. However, Civalleri et 
«Zj 18 ' 23 demonstrated the need to adjust the parameter- 
ization (from PBE-D to PBE-D*) for the application of 
the method to crystalline solids and we observe a better 
transfcrrability of the modified set for the challenging 
case of Mg(BH 4 ) 2 . 

We describe the considered library of structure types 
and the simulation settings in Sec. [TT1 present our sys- 
tematic comparison of the performance of five DFT- 
based methods for Mg(BH 4 ) 2 at T = K in Sec. Hill 
demonstrate the effect of vibrational entropy on the poly- 
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morphs' relative stability at finite temperatures in Sec. 
HVl and conclude in Sec. [Vj 



II. SIMULATION SETUP 

Library of structure types. In addition to a large pool 
of previously considered candidate structures, we include 
potentially relevant A(BC4)2 prototypes found in the In- 
organic Crystal Structure Database (ICSD^i; the full list 
of 36 structures is given in Table HI All energies are refer- 
enced to the experimental low-T a phase, which has been 
recently argued to have a higher symmetry (P6i22) 7,8 
than originally thought (P6i) 9,35 . Our simulations con- 
firm that, while the two structures are nearly degenerate 
in energy, only P6i22 is dynamically stable (P61 has mul- 
tiple imaginary phonon modes reaching 50i cm -1 ). The 
structure of the high-T /3-phase, observed experimentally 
to be stable between 453 K and 613 K 
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remains an 



open question: the first powder diffraction orthorhombic 
solution with the Fddd symmetry, based on the posi- 
tions of the B and Mg atoms, has six imaginary phonon 
modes at T reaching 65i cm -1 . Notable Mg(BH4)2 can- 
didates identified in previous DFT studies include the 
trigonal P3ml structures^ as well as the unusual low- 
density F222 phase derived from I4m 2 12 ' 14 . Figures la-c 
illustrate the diversity of the morphologies as different ar- 
rangements of the Mg 2S+ ions and BH^ - units result in 
layered, hollow-framework-like, or fairly uniform densely- 
packed structures. 

Simulation settings. The total energy calculations are 
carried out in the generalized gradient approximation 
(GGA) with the PBB-2 exchange-correlation (xc) func- 
tional and the local density functional approximation 
(LDA) 1 ^ as implemented in VASP^. We employ projec- 
tor augmented-wave pseudopotentials (PAW)- - in which 
the semicore states are treated as valence. An energy cut- 
off of 500 eV and dense Monkhorst-Pack k-point meshes 
(~ 0.03A -1 in each direction in the Brillouin zone) 38 
are applied. The vdW-DF, PBE, PBE-D, and PBE- 
D* calculations are performed with a modified version 
of PWscf 3 -^ using ultrasoft pseudopotentials. All struc- 
tures were fully relaxed with the threshold of 10 -5 [Ry] 
for energy convergence and 3T0 -1 [Ry a.u. -1 ] for resid- 
ual force; the residual stress was typically under 1 kbar. 
Benchmark tests given in Appendix A suggest that for a 
given xc functional the errors arising from other factors 
(convergence criteria, choice of pseudopotential, ect.) are 
rather small allowing us to resolve polymorphs ~ 0.005 
cV/BH 4 apart. Phonon spectra and the Gibbs energy 
corrections for selected structures are calculated with 
PBE and PBE-D* with a finite displacement method as 
implemented in PHON 40 . Forces for phonon calculations 
were obtained from VASP (PBE) or QMPOT— ^ linked 
to VASP (PBE-D*) (for more details see Sec. HVj) . 

Semiempirical dispersion contributions. The Grimme 
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FIG. 1: (Color online) a-c) Hypothetical structures of 
Mg(BH4)2 with magnesium, boron, and hydrogen atoms 
shown as large (yellow), medium (black), and small (blue) 
spheres. The isosurfaces corresponding to charge density of 
0.06 eA -3 illustrate the layered and the 3D hollow framework 
structure of P3ml and F222 phases, respectively; anisotropy 
of the charge density in I4i/acd (not shown) is less pro- 
nounced, d) Binding energy as a function of the interlayer 
spacing in P3ml-Mg(BH4)2; the layers are allowed to relax 
but the intralayer distortions are insignificant. 



corrections 24,25 are introduced as 
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where = ^/CgCg, / dmp = (1 + e-W^-D)- 1 , and 
R r is the sum of individual Rr, . 

For PBE-D calculations^ 5 - we used the following set of 
C6 coefficients [Jnm 6 mol _1 ] and van der Waals radii Rq 
[A]: for H atom C 6 = 0.14, R =1.001; for B atom C 6 
= 3.13, R Q =1.485; for Mg atom C 6 = 5.71, R =1.364. 
Dimensionless parameter in the damping function and 
global scaling parameter are d = 20 and s 6 = 0.75, re- 
spectively. For PBE-D* parameterization^ 3 - Rr, was mul- 
tiplied by 1.3 (H), or by 1.05 factor (B, Mg). 



III. PERFORMANCE OF DFT-BASED 
METHODS FOR MG(BH 4 ) 2 

The necessity to treat Mg(BH4)2 beyond the stan- 
dard DFT approximations becomes evident when the 
PBE functional shows virtually no interlayer binding for 
the exemplary layered P3ml phase (Fig. Id). The lack 
of a noticeable covalent or electrostatic interaction be- 
tween the layers is a result of the particular packing 
of the BH 4 tetrahedra that makes the interlayer inter- 
face consist of two parallel sheets of hydrogen atoms al- 
ready engaged in covalent B-H bonds. Introduction of 
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FIG. 2: (Color online) a-c) Relative Gibbs energy of selected Mg(BH4)2 polymorphs referenced to the P6i22 phase: in b) the 
relative energy is shown for all five approximations at T—0 K without zero point energy while in a,c) vibrational contributions 
to G(T) are included for PBE and PBE-D*. d) Density of Mg(BH 4 ) 2 polymorphs fully relaxed with PBE-D* and vdW-DF 
methods (p e a xp = 0.782 g cm" 3 ^. 



the vdW interactions via vdW-DF, PBE-D, and PBE-D* 
leads to interlayer cohesions of 0.123, 0.076, and 0.068 
CV/BH4, respectively. The size of the extra binding is 
substantial considering that the reported energy differ- 
ence between the experimental (a) and theoretical (F222) 
ground states within PBE is only 0.024 eV/BH^. The 
LDA, known to mimic dispersive interactions (e.g., in 
graphite^), also gives a 0.058 eV/BH4 binding, but - as 
we show below - should not be used as a substitute for 
a properly constructed nonlocal functional. 

We compare the performance of the five different meth- 
ods, PBE, LDA, PBE-D, PBE-D*, and vdW-DF, by plot- 
ting the relative stability and compound density for se- 
lected structures in Fig. 2 (the corresponding values for 
all 36 structures are listed in Table I). A number of un- 
expected features emerge as we examine the total energy 
data moving from left to right in Fig. 2b. At the PBE 
level, we find the highest stability of the least compact 
F222 phase, in agreement with Ref^. However, in the 
LDA, arguably better suited for simulating Mg(BH4)2 
based on our test in Fig. 1, there is a dramatic change 
in the ordering of the polymorphs' enthalpies with the 
ground state now being the most compact I4i/acd struc- 
ture. It is evident that, compared to PBE, the LDA 
favors higher packing: it lowers the relative enthalpy of 
the two densest P3ml and I4i/acd polymorphs by 0.134 
and 0.120 eV/BH4, respectively, and increases the rela- 
tive enthalpy of the sparsest F222 by 0.083 eV/BH 4 . 
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FIG. 3: (Color online) a — f) Nearest neighbor histograms 
in selected structures fully relaxed in PBE-D*; the stars in /) 
mark intralayer distances, g, h) Strength of pairwise Grimme 
corrections in the PBE-D and PBE-D* parameterizations. 
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FIG. 4: (Color online) Breakdown of dispersion energy con- 
tributions calculated at PBE-D* level as a function of poly- 
morph density for Mg(BH4)2 structures listed in Table [I] 



leads to a significant reordering of the relative stabilities 
obtained with PBE and promotes the I4i/acd structure. 
In fact, there is more agreement between the LDA vs 
PBE-D rather than the PBE vs PBE-D sets. The P3ml 
polymorph gets an even larger 0.217 eV/BH4 gain in rel- 
ative stability dropping below a by 0.023 eV/BH 4 , which 
may look surprising given the similar size of the interlayer 
binding for this structure in LDA and PBE-D (Fig. 1). 
To help explain this result we plot pairwise interaction 
strengths (Fig. 3g) and histograms of nearest neighbors 
(NNs) for relevant structures (Fig. 3a-f). It has been 
discusse d 23 ' 27 that introduction of the dispersive inter- 
actions 'on top' of DFT requires a smooth cut-off of the 
attractive Cij/RZ terms just below the typical vdW dis- 
tances to avoid their unphysical contributions at short 
Rij. Panels / and g in Fig. 3 show that a number of in- 
tralayer distances in P3ml (marked with stars) happen 
to fall in the critical range, around the corresponding 
ffvdw _|_ fqdw vames fitted to molecular datasets, which 
can give rise to artificially strong binding. 

The suggested adjustment of the vdW parameters 
within PBE-D* 2 ™ helps reduce the vdW overbinding for 
the dense P3ml and I4i/acd structures and for the first 
time we observe the a-phase to be the most stable poly- 
morph among the considered diverse set of candidates, 
albeit by a small margin. The classical pairwise repre- 
sentation of the vdW corrections in the Grimme approach 
(Eq. 1) makes it possible to examine the individual vdW 
contributions. According to Fig. 3h, the increase in R^ dw 
values (from 1.001, 1.485, 1.364 A to 1.301, 1.559, 1.432 
A for i — H, B, Mg, respectively) greatly diminishes the 
direct H-H interaction. However, Fig. 4 shows that the 
BH4-BH4 interaction still dominates the long-range bind- 
ing, which - not unexpectedly - scales almost linearly 
with the crystal density. The second largest contribu- 
tion, Mg-BHzi, deviates from linear dependence due to a 
variable coordination of Mg ions. The classical attractive 
terms cause up to ~ 20% volume reduction compared to 
structures optimized with PBE. 



Finally, the vdW-DF method correctly predicts a to 
be the low-T ground state and separates it from all other 
unrelated structures by a considerable 0.040 eV/BH4; the 
near degeneracy of a and P3il2 (Mn(BH4)2 prototype) 
can be traced to the close structural relationship between 
the two layered structures discussed in detail in Ref. 26 . 
The method places the hypothetical F222 structure 0.076 
cV/BH 4 above a and differentiates between the two high- 
density 2D P3ml and 3D I4i/acd structures by improv- 
ing the stability of just the latter. This suggests that the 
overestimation of the intralayer binding in the difficult 
P3ml case might still not be fully corrected within PBE- 
D*. Comparison of crystal densities for PBE-D* and 
vdW-DF in Fig. 2d shows that the latter results in less 
compact structures; the effect of the nonlocal correlations 
on the volume change is not easy to quantify because the 
method is based on the exchange and correlation func- 
tionals taken from revPBE— i and LDA, respectively. The 
calculated vdW-DF value of p — 0.758 gcm~ 3 for the a 
phase is 3% below the experimental value which is con- 
sistent with a typical 1% overestimation of bond lengths 
for GGA-based methods. 



IV. INCLUSION OF THE VIBRATIONAL 
CONTRIBUTIONS TO G(T) 

Inclusion of finite temperature contributions due to the 
vibrational entropy is orders of magnitude more expen- 
sive than the calculation of the total energy at zero tem- 
perature. We estimate the vibrational corrections using 
two DFT-based methods, PBE-D and PBE-D*, which is 
computationally feasible and allows us to show directly 
the effect of the vdW interactions. Hessian matrices were 
calculated in a numerical fashion from analitycal forces 
and atomic displacements for sufficiently large supercells 
(at least 88 atoms) 4^ According to our extensive tests 
for Mg(BH 4 ) 2 , the 0.00133 A displacements and the 10" 8 
eV SCF energy convergence criterion ensured minimal er- 
rors (~ 10 cm -1 ) for the calculated frequencies coming 
from the unharmonic effects and numerical factors, q- 
mesh for the dynamical matrices was typically two times 
denser than the fc-mesh for the energy calculations in 
each direction. The analitycal forces at the PBE level 
were found with VASP for fully relaxed polymorph struc- 
tures. The analytical forces at the PBE-D* level were 
calculated using QMPOT linked with VASP. In this case, 
since QMPOT does not currently offer stress calculations, 
the structures were optimized as follows: the unit cells 
were first fully relaxed with PWscf and then the ionic po- 
sitions were re-optimized with a combination of QMPOT 
and VASP. As described in Appendix A, our careful tests 
for selected structures showed insignificant variation in 
the equilibrium cell parameters in the two settings. Once 
the phonon densities of states n(uj) were found, the vi- 
brational contributions G v hon{T) were included via— i 

G phon (T) = k B T \n[2smh(huj /{2k B T)}n{uj)duj (2) 
Jo 
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Symm. 
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36 


Mg(AlH 4 ) 2 


P3ml (164) 


1 


0.780 


0.439 


0.875 


0.550 


0.948 


0.543 


0.735 


0.606 


0.975 


0.684 


35 


Mg(BH 4 ) 2 


P2/m (10) 


2 


0.793 


0.405 


0.915 


0.447 


1.061 


0.345 


0.807 


0.517 


1.112 


0.433 


34 


Ti T / T T \ 

Mg(BH 4 ) 2 


P3 (147) 


9 


0.845 


0.341 


0.936 


0.322 


0.985 


0.241 


0.848 


0.409 


1.037 


0.350 


33 


Zn(Re04)2 


P3 (147) 


1 


0.850 


0.340 


0.932 


0.319 


0.981 


0.238 


0.850 


0.410 


1.040 


0.350 


32 


ZrMo2 0s 


P31c (163) 


6 


0.848 


0.340 


0.934 


0.319 


0.979 


0.238 


0.850 


0.409 


1.028 


0.350 


31 


Mn(BF 4 )2 


Pnma (62) 


4 


0.847 


0.278 


0.940 


0.251 


0.993 


0.175 


0.852 


0.318 


1.046 


0.258 


30 


Mg(BH 4 ) 2 


Pnnm (58) 


2 


0.735 


0.266 


0.911 


0.270 


0.980 


0.250 


0.775 


0.320 


1.015 


0.297 


29 


Cr(AlCl 4 ) 2 


Pca2i (29) 


4 


0.809 


0.248 


0.980 


0.193 


1.051 


0.148 


0.787 


0.090 


1.085 


0.189 


28 


Mg(BH 4 ) 2 


P2i/c (14) 


2 


0.872 


0.245 


0.959 


0.173 


1.011 


0.092 


0.816 


0.186 


1.052 


0.174 


27 


Mg(BH 4 ) 2 


P2i (4) 


2 


0.859 


0.212 


1.050 


0.167 


1.117 


0.091 


0.846 


0.210 


1.149 


0.157 


26 


Ca(BH 4 ) 2 


Fddd (70) 


8 


0.903 


0.203 


0.987 


0.114 


1.043 


0.016 


0.943 


0.206 


1.090 


0.093 


25 


Mg(BH 4 ) 2 


P2/m (10) 


2 


0.963 


0.195 


1.058 


0.106 


1.101 


-0.023 


0.946 


0.193 


1.157 


0.060 


24 


Sr(AlCl 4 ) 2 


P2/c (13) 


2 


0.978 


0.194 


1.058 


0.107 


1.101 


-0.023 


0.947 


0.194 


1.158 


0.060 


23 


Mg(BH 4 ) 2 


P31m (162) 


6 


0.974 


0.194 


1.058 


0.107 


1.101 


-0.023 


0.947 


0.194 


1.153 


0.060 


22 


Mg(BH 4 ) 2 


P3ml (164) 


1 


0.968 


0.191 


1.058 


0.107 


1.101 


-0.023 


0.950 


0.194 


1.145 


0.060 


21 


ZrMo 2 8 


Pmn2i (31) 


2 


0.910 


0.182 


0.999 


0.138 


1.051 


0.080 


0.894 


0.224 


1.076 


0.143 


20 


Be(BH 4 ) 2 


I4icd (110) 


16 


0.909 


0.173 


0.878 


0.110 


0.911 


0.087 


0.677 


0.026 


0.961 


0.093 


19 


Ti T / T T \ 

Mg(BH 4 ) 2 


PT (2) 


4 


0.685 


0.172 


0.809 


0.158 


0.896 


0.149 


0.651 


0.147 


0.924 


0.175 


18 


Mg(BH 4 ) 2 


PI (1) 


2 


0.780 


0.149 


0.915 


0.111 


0.964 


0.040 


0.760 


0.086 


0.986 


0.098 


17 


Ba(BF 4 ) 2 


P2i/c (14) 


8 


0.884 


0.147 


0.980 


0.117 


1.020 


0.065 


0.874 


0.187 


1.076 


0.108 


16 


Be(A10 4 ) 2 


PT (2) 


2 


0.770 


0.147 


0.952 


0.129 


0.984 


0.103 


0.823 


0.163 


1.031 


0.082 


15 


Mg(BH 4 ) 2 


P2/c (13) 


4 


0.857 


0.140 


0.887 


0.113 


0.936 


0.101 


0.755 


0.137 


0.985 


0.120 


14 


Mg(BH 4 ) 2 


C2/c (15) 


4 


0.832 


0.113 


0.879 


0.105 


0.923 


0.095 


0.757 


0.130 


0.970 


0.114 


13 


Mg(BH 4 ) 2 


Pm (6) 


1 


0.877 


0.107 


0.950 


0.099 


0.978 


0.075 


0.841 


0.137 


1.017 


0.079 


12 


Pt(S0 4 ) 2 


P2i/c (14) 


4 


0.713 


0.106 


0.774 


0.097 


0.803 


0.088 


0.679 


0.083 


0.844 


0.102 


11 


Ca(BF 4 ) 2 


Pbca (61) 


8 


0.841 


0.091 


0.912 


0.085 


0.958 


0.063 


0.766 


0.105 


1.008 


0.072 


10 


Mg(BH 4 ) 2 


I4m2 (119) 


4 


0.529 


0.077 


0.553 


0.041 


0.553 


0.048 


0.533 


-0.019 


0.585 


0.041 


9 


Mg(BH 4 ) 2 


F222 (22) 


8 


0.522 


0.076 


0.552 


0.040 


0.552 


0.047 


0.532 


-0.020 


0.591 


0.063 


8 


Cd(AlCl 4 ) 2 


Pc (7) 


2 


0.836 


0.062 


0.936 


0.066 


0.967 


0.052 


0.801 


0.092 


0.998 


0.054 


P 


Mg(BH 4 ) 2 


Fddd (70) 


64 


0.775 


0.062 


0.790 


0.056 


0.819 


0.044 


0.727 


0.037 


0.865 


0.049 


6 


Mg(BH 4 ) 2 


Pmc2i (26) 


2 


0.835 


0.058 


0.934 


0.066 


0.966 


0.052 


0.790 


0.090 


1.002 


0.054 


5 


Mg(BH 4 ) 2 


I4i/amd (141) 


8 


0.939 


0.047 


1.030 


0.018 


1.044 


-0.007 


0.974 


0.103 


1.082 


-0.007 


4 


Sr(AlCl 4 ) 2 


I4i/acd (142) 


4 


0.972 


0.043 


1.049 


0.002 


1.078 


-0.027 


0.986 


0.098 


1.102 


-0.022 


3 


Sr(AlCl 4 ) 2 


Pbca (61) 


8 


0.944 


0.040 


1.047 


0.004 


1.074 


-0.023 


0.979 


0.098 


1.100 


-0.019 


2 


Mn(BH 4 ) 2 


P3il2 (151) 


9 


0.780 


0.000 


0.828 


0.007 


0.837 


0.003 


0.791 


0.010 


0.878 


0.004 


a 


Mg(BH 4 ) 2 


P6i22 (178) 


30 


0.758 





0.808 





0.814 





0.770 





0.860 






TABLE I: Relative stability AE [eV/BH 4 ] and density g [g cm -1 ] (T = K, no zero point energy corrections) of magnesium 
borohydride polymorphs referenced to the experimentally observed a phase. Z stands for the number of Mg(BH 4 ) 2 formulas 
per unit cell. In the third column space groups information (international symbols and numbers) is reported. 



Calculated relative Gibbs energies G(T) for selected 
low-energy structures are shown in Fig. 2a and 2c. We 
find that the zero point energy alone can change the rel- 
ative stability by over 0.04 eV/BH 4 . Addition of dis- 
persive interactions leads to about 0.01 eV/BELi changes 
in the relative stability at T = 800 K and can affect 



the ordering of the polymorphs as happens for the P3il2 
structure. The related a and P3il2 phases remain nearly 
degenerate in the whole T range and cannot be unam- 
biguously resolved within our vdW-DF or PBE-D* sim- 
ulations. The search for the high-T ground state may be 
simplified by our observation of little variation in the en- 
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thalpy difference between a and /3 for all the five methods 
in Fig. 2: due to the close structural relationship between 
the two phases, shown in Fig. 3a and 3c, there must be 
a cancellation of errors. Hence, despite the remaining 
noticeable difference between the PBE-D* and vdW-DF 
relative stability for some polymorphs the former could 
be an appropriate choice for identification of viable high- 
T candidates and estimation of the phonon contributions 
to G(T). 

V. CONCLUSIONS 

The presented results illustrate that identification of 
ground states depends not only on the exhaustive sam- 
pling of possible structures and compositions but also 
on the accuracy of the chosen simulation method. We 
have demonstrated that inclusion of the dispersive inter- 
actions via nonlocal correlation functional within vdW- 
DF leads to a good agreement with experiment for 
the Mg(BH4)2 low-T ground state. The related well- 
characterized Ca(Bfi4)2 system should be an interesting 
test case for the vdW-DF and PBE-D* methods since 
standard DFT approximations seem to already reproduce 
the experiment^. Proper treatment of the nonlocal cor- 
relation effects is expected to be most critical in higher- 
valent metal borohydrides that can form weakly interact- 
ing molecular-type complexes^. Until such large crys- 
talline systems can be handled with accurate quantum 
Monte Carlo-level or quantum chemistry methods, the 
continued improvement of DFT-based methods^ may 
provide an attractive alternative for description of the 
dispersive interactions therein. In any case, it appears to 
be good practice to systematically test proposed ground 
states with a range of DFT flavours as this can help spot 
and avoid potential artifacts of DFT approximations. 
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Appendix A: Benchmark tests 

To estimate how technical factors, such us a choice of 
a pseudopotential, may influence the accuracy of total 
enthalpy calculations we performed an optimization of a 
set of exemplary structures at the PBE-D* level using the 
projector augmented- wave pseudopotentials (PAW)2£, as 
implemented in VASP, as well as the ultrasoft pseudopo- 
tentials (USPP)^, as implemented in PWscf (see Table 
HE])- In the first case, we used QMPOT software linked to 
VASP to calculate dispersion contributions to energy and 
forces. A crystal cell has been optimized in a series of in- 
dependent runs for a set of cell parameters. For all exem- 
plary structures the relative difference in cell parameters 
optimized using PAW pseudopotential (VASP+QMPOT) 
and USPP pseudopotential (PWscf software) was below 
0.5%. The difference between relative energies calcu- 
lated with PAW and USPP was found to be below 0.0035 
eV/BtLi. We estimate that for a given xc functional the 
errors arising from other factors (convergence criteria, 
choice of pseudopotential, ect.) are rather small allowing 
us to resolve polymorphs ~ 0.005 eV/BH 4 apart. 





a 


2 


3 


4 


5 


6 


8 


10 


11 


12 


22 


PAW 
USPP 
difference 


0.0000 
0.0000 
0.0000 


0.0065 
0.0065 
0.0000 


0.0029 
0.0041 
-0.0012 


0.0022 
0.0022 
0.0000 


0.0188 
0.0179 
0.0009 


0.0670 
0.0660 
0.0010 


0.0671 
0.0658 
0.0013 


0.0373 
0.0406 
-0.0033 


0.0848 
0.0845 
0.0003 


0.0965 
0.0968 
0.0003 


0.1091 
0.1066 
0.0025 



TABLE II: Relative enthalpy [eV/BHU] of selected polymorphs calculated at the PBE-D* level with different pseudopotentials. 
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